library(tidyverse)
library(plotly)
library(sf)
library(mapview)
library(tigris)
library(censusapi)
library(readxl)
library(leaflet)


options(
  tigris_class = "sf"
)

Sys.setenv(CENSUS_KEY="10dcd73d7c043e91bac9fb8d3989cbff54b08790")

Mapping populations of individuals 65 and older relative to hospital bed counts

Here I examine the total number of hospital beds in facilities with the designation “Acute Critical Care Hospital” from the California Department of Public Health data, per 1000 members of the population of age 65 and older in each county in California.

First we can look generally at the percent of the population that is of age 65 and older by county in California.

# If not already done, obtain the census data
# acs_vars <-
#   listCensusMetadata(
#     name = "2018/acs/acs5",
#     type = "variables"
#   )
# # Save an .rds version on my computer
# setwd("~/Documents/2020 Spring Quarter/CEE 218Z")
# saveRDS(acs_vars, file = "censusData2018_acs_acs5.rds")
# setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")

# obtain the saved census data 
setwd("~/Documents/2020 Spring Quarter/CEE 218Z")
acs_vars = readRDS("censusData2018_acs_acs5.rds")
setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")

# get sex by age for CA counties, modified from the rBasics.Rmd code
CA_county_sexbyage <-
  getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "county:*", 
    regionin = "state:06",
    vars = "group(B01001)"
  ) %>%
  mutate(
    countyval =
      paste0(state,county)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
  dplyr::select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
  gather(
    key = "variable",
    value = "estimate",
    -countyval
  ) %>%
  mutate(
    label = acs_vars$label[match(variable,acs_vars$name)]
  ) %>% 
  select(-variable) %>% 
  separate(
    label,
    into = c(NA,NA,"sex","age"),
    sep = "!!"
  )


# Now find the population that is age 65 and up - also modified from rbasics.Rmd
CA_county_elderly <- 
  CA_county_sexbyage %>% 
  filter(!is.na(age)) %>% 
  mutate(
    elderly = 
      ifelse(
        age %in% c(
          "65 and 66 years",
          "67 to 69 years",
          "70 to 74 years",
          "75 to 79 years",
          "80 to 84 years",
          "85 years and over"
        ),
        estimate,
        NA
      )
  ) %>% 
  group_by(countyval) %>% 
  summarize(
    elderly = sum(elderly, na.rm = T),
    total_population = sum(estimate, na.rm = T)
  ) %>% 
  mutate(
    percent_elderly = elderly/total_population
  )


# get the counties as sf objects in CA
CA_counties <- counties("CA", cb = F, progress_bar=F)

# map percent elderly
CA_county_elderly %>% 
  left_join(CA_counties %>% dplyr::select(GEOID), by = c("countyval" = "GEOID")) %>% 
  st_as_sf() %>% 
  filter(!is.na(percent_elderly)) %>% 
  mapview(zcol = "percent_elderly", layer.name = 'Percent ages 65+')

Now let’s combine this data with the hospital data and map the total number of beds in Acute Critical Care Hospitals per 1000 members of the population ages 65 and above.

# get data on hospital beds
setwd("/Users/simonespeizer/pCloud Drive/Shared/SFBI/Data Library/CHHS")
hospBeds <- read_excel("healthcare_facility_beds.xlsx")
setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")

# filter by general acute care hospitals
gen_acute_beds <- hospBeds %>% filter(FAC_FDR == 'GENERAL ACUTE CARE HOSPITAL') %>%
  group_by(FACNAME, COUNTY_NAME, BED_CAPACITY_TYPE) %>% summarize(gen_acute_beds_hosp_type = sum(BED_CAPACITY, na.rm = T)) 

# number of beds at each hospital 
gen_acute_beds_by_hosp <- gen_acute_beds %>% ungroup() %>% group_by(FACNAME, COUNTY_NAME) %>% 
  summarize(total_acute_beds_per_hosp = sum(gen_acute_beds_hosp_type, na.rm = T))

# number of beds in each county
gen_acute_beds_county <- gen_acute_beds_by_hosp %>% ungroup() %>% 
  group_by(COUNTY_NAME) %>%
  summarize(totalHosp = n(), total_acute_hosp_beds = sum(total_acute_beds_per_hosp, na.rm = T)) 

# add county names in upper case to be able to join
CA_counties <- CA_counties %>% mutate(uppercasename = toupper(NAME))

# join county data and beds data
gen_acute_beds_county_w_geom <- CA_counties %>% left_join(gen_acute_beds_county, by = c('uppercasename' = 'COUNTY_NAME')) %>% mutate_all(~(ifelse(is.na(.), 0, .)))

# join county data and elderly data and get beds per population and per elderly population
gen_acute_beds_county_final <- gen_acute_beds_county_w_geom %>%
  left_join(CA_county_elderly, by = c("GEOID" = "countyval")) %>%
  mutate(beds_per_1000_elderly = 1000*total_acute_hosp_beds / elderly, beds_per_1000 = 1000*total_acute_hosp_beds / total_population)

# map beds per elderly population
gen_acute_beds_county_final %>% 
  dplyr::select(beds_per_1000_elderly) %>% 
  mapview(layer.name = "Acute critical care hospital beds per 1000 residents ages 65+")

Let’s also look at the number of acute critical care hospital beds per 1000 members of the total population, just for comparison.

gen_acute_beds_county_final %>% 
  dplyr::select(beds_per_1000) %>% 
  mapview(layer.name = "Acute critical care hospital beds per 1000 residents")

These numbers include some beds that may not realistically be available readily or useful for COVID-19 patients, such as intensive care newborn beds. Let’s also consider the number of beds specifically for intensive care, coronary care, acute respiratory care, intermediate care, rehabilitation, and unspecified general acute care, assuming that these are the ones most likely to be available for adult COVID-19 patients.

# define the bed categories to be considered
selected_bed_cats <- c('ACUTE RESPIRATORY CARE', 'CORONARY CARE', 'INTENSIVE CARE', 'INTERMEDIATE CARE', 'REHABILITATION', 'UNSPECIFIED GENERAL ACUTE CARE')

# get the selected bed types
selected_gen_acute_beds <- gen_acute_beds %>% 
  filter(BED_CAPACITY_TYPE %in% selected_bed_cats)

# number of beds at each hospital 
selected_gen_acute_beds_by_hosp <- selected_gen_acute_beds %>% ungroup() %>% group_by(FACNAME, COUNTY_NAME) %>% 
  summarize(total_acute_beds_per_hosp = sum(gen_acute_beds_hosp_type, na.rm = T))

# number of beds in each county
selected_gen_acute_beds_county <- selected_gen_acute_beds_by_hosp %>% ungroup() %>% 
  group_by(COUNTY_NAME) %>%
  summarize(totalHosp = n(), total_acute_hosp_beds = sum(total_acute_beds_per_hosp, na.rm = T)) 

# join county data and beds data
selected_gen_acute_beds_county_w_geom <- CA_counties %>% left_join(selected_gen_acute_beds_county, by = c('uppercasename' = 'COUNTY_NAME')) %>% mutate_all(~(ifelse(is.na(.), 0, .)))

# join county data and elderly data and get beds per population and per elderly population
selected_gen_acute_beds_county_final <- selected_gen_acute_beds_county_w_geom %>%
  left_join(CA_county_elderly, by = c("GEOID" = "countyval")) %>%
  mutate(beds_per_1000_elderly = 1000*total_acute_hosp_beds / elderly, beds_per_1000 = 1000*total_acute_hosp_beds / total_population)

# map beds per elderly population
selected_gen_acute_beds_county_final %>% 
  dplyr::select(beds_per_1000_elderly) %>% 
  mapview(layer.name = "Selected types of acute critical care hospital beds per 1000 residents ages 65+")
# map beds per total population
selected_gen_acute_beds_county_final %>% 
  dplyr::select(beds_per_1000) %>% 
  mapview(layer.name = "Selected types of acute critical care hospital beds per 1000 residents")